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Abstract 

Starting from Enskog equation of hard spheres of mass m and diameter 
D under the gravity g, we first derive the exact equation of motion for the 
equilibrium density profile at a temperature T and examine its solutions via 
the gradient expansion. The solutions exist only when (3fj, < /i ~ 21.756 in 
2 dimensions and // D ~ 15.299 in 3 dimensions, where \x is the dimensionless 
initial layer thickness and (5 = mgD/T. When this inequality breaks down, a 
fraction of particles condense from the bottom up to the Fermi surface. 

PACS numbers: 05.20-y, 81.35+k, 05.20. Dd, 05.70.Fh 

Granular materials are basically a collection of hard spheres that interact with each other 
via hard sphere potential [1]. For this reason, many of the properties of excited granular 
materials may be understood from the atomistic view of molecular gases, in particular 
from the view point of kinetic theory [2] . There are, however, several distinctions between 
molecular gases and granular materials: First, granular materials are macroscopic particles 
with finite diameter, and thus they cannot be compressed indefinitely. Second, the gravity 
plays an important role in the collective response of granular materials to external stimuli, 
largely because of the ordering of grains induced by the gravity. For example, one of the 
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notable characteristics of the excited granular materials in a confined system under gravity 
is the appearance of a thin boundary layer near the surface that separates a fluidized region 
from a rigid solid region. This has been known for long time as is evidently seen in shearing 
experiments [3], avalanches [4], and grains subjected to weak excitations [5,6]. In this limit, 
those grains in a solid region are effectively frozen, and thus do not participate in dynamical 
processes. Now, the kinetic theory relies on two particle collision dynamics, and has been 
applied to systems where all the granular particles are in motion colliding with each other. 
But if a portion of grains are frozen and remain largely motionless, the kinetic theory or in 
general the continuum theory may pose some problems. Consider for example a system of 
strongly excited granular particles under gravity, where all the grains undergo collisions and 
thus the kinetic theory is valid. If we decrease the strength of excitation, then the particles 
at the bottom will freeze themselves, and the boundary layer will develop at the top. The 
question we address in this Letter is: How does the kinetic theory describe such process? 

In a recent paper [6] , it has been demonstrated that the granular statistics in the presence 
of gravity does not follow the usual Boltzmann statistics as in molecular gases, where all the 
particles are dynamically active, but a new Fermi statistics, where most of the particles are 
effectively frozen and only a portion of particles near the surface participate in the dynamical 
process. This is due to the excluded volume effect and the ordering of potential energy by 
gravity, and the mechanism associated with this Fermi statistics is similar to that of the 
Fermi gas in a metal. The existence of a thin boundary layer in granular materials should 
be viewed from such perspective. 

Our specific objective of this Letter is to use the kinetic theory, in particular the Enskog 
equation of hard spheres of mass m and diameter D, to explore whether or not the kinetic 
theory can describe the cross over from Boltzmann to Fermi statistics and if so, under what 
conditions it occurs. Our particularly interesting discovery is that the prediction of the 
Enskog equation is only valid when Pfj, < where /i is the dimensionless initial layer 
thickness of the granules(or the Fermi energy), (3 = mgD/T with T the temperature, and 
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the critical value, // , is determined to be /j, = 21.756 in 2d and jjt, = 15.299 in 3d. When 
this inequality is violated, Enskog equation does not conserve the particles, and the missing 
particles condense from the bottom up to the Fermi surface. This way, the hard sphere 
Enskog gas appears to contain the essence of Fermi statistics and Bose condensation. 

The starting point of our investigation is the Enskog's kinetic equation for elastic hard 
spheres [7] of mass m and diameter D in the presence of gravity: 



- + v.Vf-m 9 — = J E (1) 



where the Enskog's collision operator is given by: 



Je = D 2 J d 3 Vl j + d 2 ee.g[f(v, V)f(r+De, v' 1 ) X (r+ ±De) -f(v, v)/(r-£>e, Vl ) X (r- ^De) 

(2) 

Here, (v, vi) and (v', v[) are the velocities of two colliding particles before and after the 
collision, and r and r' are the positions of the two particles when they are in contact, e is 
the unit vector in the direction of r-r', V is the relative particle flux given by V = v — v x 
and V = v' — Vi' [2]. The + sign means that the integration should be carried out with 
the restriction that < e • V. From the geometry it is easy to show V = V — 2ee • V and 
|V|' = |V|. The correlation function, x, is defined as F 2 (r, v; r 1; v x ) = i*i(r, v)F 1 (r, Vi)x, 
where F 2 and F 1 are the two and one particle distribution function respectively. Molecular 
chaos assumption in the usual Boltzmann collision operator sets x — 1- 111 the case of 
dense gases, Enskog [2,7] assumed that x might be given by the equilibrium two point 
correlation function estimated at the contact point. At equilibrium, we expect that the 
factorization of the space and velocity in the distribution function is valid. Thus, we set: 
f(r,v;t = oo) = G(r)0(v). This is equivalent to the separation of the configurational 
statistics from the kinetics. In equilibrium statistics of elastic granular materials, 0(v) is 
expected to be Gaussian, and the energy conservation requires that 0(v' 1 )0(v / ) = 0(vi)0(v). 
We note that 0(v) is normalized to one: / (i 3 v^(v) = 1. In the steady state, the left hand 
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said of eq.(l) then becomes: 



^(v)[^ + ^G(,)] (3) 

In order the compute the right hand side, we first notice that the collision operator 
J E (eq.(2)) is given by: 



J £ = J D0(v)G(r)[|rf 2 v0(v 1 )] J d6e\x(T~De)G(T-De)- X (T+^De)G(T+De)]e»V (4) 

We now first compute the right hand side of eq.(l) under the factorization assumption. In 
this case, a care must be taken because the integral should be performed with the constraint 
e • V > 0. However, if we change e to — e of the second term in the Enskog operator 
Je(eq.(4)), we find it becomes: 



-/ d9 X (r + ^De)G(r + £>e)e • V = / d9e • V X (r - ^De)G(r - De) (5) 

Je»V>0 2 Je»V<0 2 

Hence, we can now remove the restriction e • V > in eq.(4) and integrate over the 
whole space. After some algebra, we obtain: 

J E = D<j>(v)G(r)[j rf 2 Vl 0( Vl )]v . I (6) 

where 

I = f d9e X (r - i J De)G(r - De) 

J allspace 2, 

= \f d9ee[ X (r - i,De)G(r - De) - X (r + ]-De)G(r + De)] (7) 

Z J allspace Z Z 

In obtaining (7), we utilized the fact that / d 2 vi4>(v 1 )v 1 = by symmetry. Our next 
step is to compute I. To this end, we first rewrite I: 

I = J d9e X (r - l -De)G(v - l -De) = - J d9e X (r + l -De)G{v + l -De) 

= \ J d9e[ X (r - l -De)G(r - De) - X (r + ^£>e)G(r + ^De)] (8) 
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Notice that I x (and I y ) will vanish by symmetry and only the vertical component, I z , 
survives. Using e = (sin 9, cos 9) in 2d, we find: 



1 r 27T 1 1 

I z = - dOcosOhrtz - -Dcos9)G(z - DcosO) - x(z + -Dcos9)G(z + DcosO)} (9) 

2 Jo 2 2 

At equilibrium, 0(v; T) = 0(mv 2 /2T). We now put this functional form along with 
G(z;T) = G{mgz/T) to the right hand side of (1) and cancel v z (f>(v) term. Next, in a free 
volume theory, particles are confined in a cage. Hence, if we use a simple cubic lattice as a 
basic lattice, the closed packed volume fraction p c = N/V = N/D 2 N = 1/D 2 . If we define 
the dimensionless density 4>(z) = G(z)/p c = D 2 G(z) or 0(C,/3) = D 2 G(z) with ( = z/D, we 
then obtain the exact dimensionless equation of motion for 



^ + /30(C) = 0(C)/ f (C) (io) 

where 



'c(C) = \ r d6 ^ - l cos6 )^ - cosd ) - X(C + \cos9m + cos6)] (11) 
Z Jo z z 

For 3d,the corresponding equation for the density <p(() = D 3 G(z) is given by: 

^ + /30 = 0/c(C) (12) 

with 

U() = 7T r d6sin6cos6[x(C ~ cos6/2)(f)(( - cosff) - X (( + cos6/2)(f)(( + cos6)} (13) 
Jo 

Several forms for the equilibrium correlation function x have been proposed, but we use 
the following widely used forms: For 2d, we use the form proposed by Ree and Hoover [8]: 
x{4>) = (1 — «i0 + «20 2 )/((l — «0) 2 , while for 3d, we use the form suggested by Carnahan 
and Starling [9]: X (4>) = i 1 ~ tt0/12)/(1 - 7T0/6) 3 

Since the total number of particles, N, remain fixed, the following normalization condi- 
tion should be satisfied for both 2d and 3d. 



poo 

/ dC<P(C P) = fi (14) 

Jo 

where \i = N/Q x (or ji = N/Q x Q y in 3d) is the Fermi energy [6] and Cl x ,£l y are the 
degeneracies along the x and y axis. We now perform the gradient expansion of (11) and 
(13) and retain only the terms to first order in dx/d(. We find: 



* +w= _^ +2x |, m (15a) 



The solutions are readily obtained. For 2d, we find: 

—(3(C ~ P) — l n 4> + c i0 + C2log(l — atf) + c 3 / (1 — atf) + c 4 / (1 — atf) 2 (16a) 

Pfi = ln<f) + ci0 o + c 2 ln(l - ct(j) ) + c 3 / (1 - a ) + c 4 / (1 - a0 o ) 2 (166) 

/3fi = (f) + ci0 2 /2 + c 2 (0 o + - atj) )/a) + c 3 /n(l - a0 o )/a 

-(c 4 /a)[l/(l - a0 o ) - 1] + c 3 o /(l - a0 o ) + c 4 o /(l - a0 o ) 2 (16c) 

where O is the density at ( — 0, and ci = 2a 2 /a 2 | ps 0.0855, c 2 = — f (c*i — 2a 2 /a)/a: 2 
0.710, c 3 = -c 2 ,c 4 = f (1 - ai/a + a 2 /a 2 /)a ps 1.278. For 3d, we find: 

-ji) = Intp - 1/(1 - atf) 2 + 2/(1 - atpf (17a) 

#i = ln{<j> ) - 1/(1 - a0 o ) 2 + 2/(1 - a<p o f (176) 



PH = $o- Z — + -7- -TTg 17C 

1 - «0 O (1 - a0 o y 
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where a = n/6. For given values of (5 and /i, (p Q = <p(( = 0) will be determined by eq.(16c) 
and (17c). However, since the right hand sides are monotonically increasing functions for 
O , fin must have the upper bound fi ; namely, /i D = 21.756 in 2d and /i = 15.299 in 3d, 
which are the values obtained by setting 4> Q = 1 in the right hand side of (16c) and (17c). 
Considering the fact that both the temperature T and the Fermi energy \i are arbitrary 
control parameters, the existence of such bounds is a puzzle: if ftp is less than /i , then the 
density profiles given by Eq.(16a) and (17a) are well determined, but if is greater than 
then O must be one, and the particle conservation breaks down, namely 



The central question is: where does the rest of the particles go? In order to gain some 
insight into this question, consider first the case of point particles under gravity, in which 
case the density profile is given by: p(() = p(0)exp(—mg(/T). If we put more particles into 
the system, we simply need to increase p(0) because the point particles can be compressed 
indefinitely, and the profile simply shifts to the right. We now replace these point particles 
with hard spheres, which cannot be compressed indefinitely. Suppose we start from a high 
temperature where all the particles are active. We then slowly decrease the temperature to 
suppress the thermal motion. At a certain temperature, the freezing of the particles will 
occur from the bottom [11], which will then spread out as the temperature is lowered down 
further, until at T=0 all the particles are frozen. Note that the frozen particles in the closed 
packed region behave like a solid. Such observation helps us to resolve the puzzle associated 
with the disappearance of particles, which must then condense from the bottom up to the 
lower part of the fluidized layer. We term this surface, which separates the frozen or a 
closed packed region with = 1 from the fluidized region with < 1, the Fermi surface. 
The location of the Fermi surface, (p, is determined by the amount of the missing particles, 
namely, (f = /-t — fi Q / (3. For nonzero (f, we must put the missing particles below the Fermi 
surface and shift the bottom layer from ( = to (p- Such modified profile for j3 — 1 and 




(18) 
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fj, = 100 in 2d is shown in Fig.l. In comparison, in Fig.l is also shown the density profile for 
P — 10 and pi = 100 for which ft/i < p: . We conclude that while the density profile obtained 
this way is not exactly the same as the Fermi profile, the essence of the Fermi statistics, 
namely the effect of excluded volume interactions, sets in when f3p > /i . The condensation 
of particles from the bottom in one dimensional vibrating bed [12] and the clustering of 
particles near the bottom wall in two dimensional experiments [13] appear to be a strong 
confirmation of this scenario, which seems to be reminiscent of the Bose condensation of 
particles into the ground state. 

We conclude with a few remarks. First, it remains to be seen whether the feature observed 
in this paper for hard spheres persists in the presence of dissipation, namely when particles 
collide inelastically. In this case, there are some evidences that the velocity distribution 
function is not Gaussian for strong dissipation, and the factorization assumption may not 
be valid. The future studies must focus on the accurate determination of this velocity profile, 
based on which an extension of the present analysis must be carried out. Second, we may 
define the freezing temperature T c as the point where the particle conservation breaks down, 
namely 

T c = mgDp/p (19) 

which may be tested experimentally presumably by using the relation between T and the 
vibration strength T [6], or by Molecular Dynamics simulations. Finally, we point out that 
eqs. (15a) and (15b) can also be obtained from the force balance eq. for the pressure 
P:dP/dz = —pg with P obtained by the virial expansion of hard sphere gas in the absence 
of gravity [10]: P/Tp = 1 + A d pD d X {D) /2d, where A d = 4tt for 3d and A d = 2tt for 2d. 

The author is particularly grateful to A. J. McLennan for numerous discussions on kinetic 
theory and for several incisive suggestions over the course of this work. The author also 
wishes to thank H. Hayakawa for helpful discussions, S. Luding for ref. [12] and Joe Both 
and Paul Quinn for checking some of the algebras. 
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Figure Caption 

Fig.l. The crossover from Boltzmann to Fermi statistics as the temperature is lowered. 
For a given Fermi energy, /i = 100, the density profiles are shown as a function of dimen- 
sionless height ( for (3 = 1/10 and (3 = l(dotted line) for the two dimensional Enskog gas. 
For (5 — 1, grains freeze from the bottom up to the Fermi surface, and only those grains 
near the Fermi surface participate in the dynamical process. 
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